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ABSTRACT 


Previous experimental investigations have shown that the characteristics of flow 
about a circular cylinder immersed in a time-dependent flow exhibit cycle-to-cycle 
variations. These variations have been attributed to the variations in the spanwise 
coherence, aspect ratio, nonuniformity of the flow, and random disturbances in the 
ambient flow. A theoretical investigation was undertaken to examine the stability of 
the flow characteristics in terms of the initial state of the vortices. An idealized model 
has been devised and the position of the vortex was varied systematically. The results 
have shown that finite-precision information about the characteristics of the flow does 
not lead to finite-precision information at a later stage. In fact the advection of the 
vortices can give rise to chaotic behavior in the calculated lift and drag forces and in 
the velocity field. It is concluded that the cycle-to-cycle variations are not entirely due 
to lack of spanwise coherence and that they are mostly a consequence of the chaotic 


motion which can result from the advection of the vortices in a time-dependent flow. 
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I. INTRODUCTION 


In the past ten years there has been much effort to analyze the separated time 
dependent flow about bluff bodies in general and circular cylinders in particular. 
Shortfalls in the Morison Equation to adequately predict in-line forces as well as the 
advent of sophisticated computational and experimental techniques have stimulated a 
great deal of research in this area. The accumulation of this numerical, experimenta! 
and visual data has generated the need for a deeper understanding of the motion of 
vortices in time dependent flow. Specifically in oscillating flow the stability of the 
vortices generated is Of particular importance. The significance of the stability of the 
vortices is a result of its possible exposition of documented cycle-to-cycle variations in 
measured forces on the cylinder. 

In this application, stability is not a characteristic of the vortex life or demise but 
of the vortex motion or path. A system is referred to as stable in any particular state 
of equilibrium if, after being given a finite perturbation, it tends to return to the state 
of equilibrium existing before it was disturbed, On the contrary an unstable system 
subjected to a finite perturbation catastrophically fails to return to its initial state of 
equilibrium. Many real systems in nature do not fall into either of these two categories 
though. In these systems there are certain times when finite disturbances do not lead 
to finite results nor catastrophic or runaway reactions. These complicated systems are 
described as neither stable nor unstable but as systems susceptible to chaos. 

Whether or not the stability of the vortex motion in oscillating flows can be 
categorized as chaotic requires not only an analysis of the vortex motion but also a 
clear definition and understanding of chaos. Unfortunately the pursuit of this 
understanding remains the subject of ongoing research and a concise definition is yet 
elusive. Strides have been made however in identifying chaotic behavior. Svstems 
studied that demonstrate chaotic behavior could best be described as bifurcating 
systems that operate along a razor’s edge or precipice. The slightest disturbance will 
push the system one way or the other and will lead to decidedly different (but not 
unstable) results. Additionally, the systems studied reveal that the state of chaotic 
behavior is not all encompassing but exists more in regions of susceptibility to chaotic 


behavior preceded and followed by regions of stable behavior. 


1] 


The application of this analysis of chaotic behavior to vortex motion in 
oscillating flow is particularly useful due to the complexity of the factors driving the 
vortex trajectory. The coupling of relative position to other vortices (real or image) in 
a time dependent flow is both intricate and inseparable. Additionally, the flow field 
about a bluff body near the surface, combined with the problems of multiple vortices, 
viscous dissipation, wall reflection and boundary layer and limited coherence length, 1s 
extremely intricate. Consequently, the motion of the vortex at times becomes subject 
to finite perturbations without finite results. The specific coupling that makes the 
system susceptible to chaotic behavior gives rise to the regions of chaos preceded and 
followed by regions of stability. 

In order to understand these complicated causes and effects first an 
understanding of the motion of a single potential vortex in oscillating flow about a 
bluff body is required. Therein lies the impetus of this study. In an ideal two 
dimensional potential flow, wall reflection. wall boundary layers and coherence length 
are elimunated. Additional vortices can be added to study mutual induction effects. 
Further, advanced models of the vortex (i.e. Rankine and Gaussian time dependent 
models) can clarify the broad effects of viscous dissipation. 

The scope of this study includes the potential flow of one and two vortices placed 
in the vicinity of a circular cylinder in oscillating flow. Emphasis is placed on 


identifying the onset of chaotic behavior and its effects on the flow characteristics. 
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Hl. BACKGROUND 


The concept of chaos in dynamic systems is by no means new. However, the 
detailed investigation as to the onset and mechanisms of chaos is still comparatively 
young. The primary difficulty in the investigation of chaos is the fact that by definition 
chaotic behavior is non-integrable. Because of this the first attempts to define the 
concept of chaos were undertaken by mathematicians. In 1984 Holmes [Ref. 1] 
investigated chaotic motion in forced oscillations. Here the onset of chaos was linked 
to periodicity or aperiodic behavior of the function through scrutiny of the harmonics. 
Although this is a credible definition its application appears limited to the idealized 
mathematical system from which it was developed. In 1980 Aref [Ref. 2] expanded on 
this definition and the use of Poincare’ maps by illustrating chaotic behavior in paths 
of finite, ideal, and controllable vortices. His motivation for the investigation was the 
onset and mechanisms of chaos. Therefore, while the idealized system that he designed 
was useful for demonstrating chaotic behavior, its application to real world dynamics 
was severly limited. 

Other efforts were made to trv to link the study of chaos to systems more closely 
resembling those found in nature. In 1984 Hardin & Mason (Ref. 3] investigated the 
two dimensional svstem of axial vortices in a pipe. This was the first attempt to 
examine a system that resembled real phenomena and that also underwent chaotic 
behavior. All previous studies examined chaotic behavior by generating systems 
susceptible to chaos and then extrapolating the results to possible real world problems. 

Finally, in 1986 Dowell & Pezeshki [Ref. 4] investigated the one dimensional 
Duffings Equation and the onset of chaos. In its relation to real systems Duffings 
Equations could be used to describe the snap buckling of a beam and other 
phenomena. In the beam the onset of chaos is easily observed but the integration of 
the equation to predict or study the mechanisms of chaos is somewhat abstract. Of 
particular interest in this study was the examination of the onset of chaos by Phase 
Plane Plots (displacement - velocity plots). 

By reasons of their treatment of fluid dynamics in general and vortex motion in 
particular, Aref's chaotic advection and Hardin & Mason’s vortices in a pipe have 


direct application to this investigation and will be discussed in greater detail. 
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A. CHAOTIC ADVECTION | 

Aref [Ref. 2] devised a numerical model to represent the motion of an 
incompressible inviscid fluid within a circular bounding contour. A stirring agitator 
was modelled as a potential point vortex. The agitator(s) provided the potential flow 
within the boundary. A global Eulerian view of an array of markers was traced to 
examine and verify the poor efficiency of the laminar mixing. 

In a follow-up study, he examined the case of two fixed agitators (point vortices) 
in the fluid that generated “piecewise-constant stirrer motion.” The vorticity of the 
agitators was controllable (on or off) and switched back and forth at prescribed times 
to generate unsteady flow. Again a marker array of particles was traced. The mixing 
was clearly superior and Aref classified the motion as chaotic. The parameters that 
dictate the switching of the vortices and therefore the chaotic regions were cited as 
mechanisms for chaos. 

Extrapolations to real phenomena include enhanced mixing chambers with 


controllable agitators and similarities to gravitational systems. 


B. VORTICES IN A PIPE 

The motion of axial vortices was investigated by Hardin and Mason [Ref. 3} 
because of its widespread appearance in physical situations. These vortices occur in 
virtually every flow that changes axial directions (1.e. rivers, pipes,.lung bronchial tubes, 
etc.). In these physical systems pairs of counterrotating vortices are produced. This 
study investigated the effects of one and two pairs of equal strength counterrotating 
two dimensional potential vortices enclosed in a circular boundary. Vortex advection 
by induced velocities (real or image) was determined and displayed bv means of 
Eulerian paths of the vortices. 

Controlling factors on this investigation were initial locations of the vortices 
within the boundary. For the case of single pair of vortices the paths were found to be 
periodic or stable. In the case of two-vortex pairs (four vortices within the boundary), 
the periodicity of the vortex path and flow was a function of where the vortices were 
placed initially. This supports the concept of chaotic motions occurring within regions 
of susceptibility. It was found that within these regions of susceptibility that small 
perturbations (0.001 radii of the enclosing boundary) were sufficient to drive the vortex 
path aperiodic and demonstrate chaotic behavior. Hardin & Mason [Ref. 3] do note 


that experimental evidence indicates that bifurcating -or branching flows are 
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quantitatively more stable in real viscous flows than these numerical results suggest. 
This is attributed to the complexity of this secondary flow as well as viscous 
dissipation, vortex decay and coring influences. It does not, however, diminish the 


validity of the study or the qualitative results of the susceptibility of fluid flows to the 
onset of chaos. 
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lit, ANALYSIS 


A. INTRODUCTION 

In 1976 Sarpkaya [Ref. 5] presented extensive data on oscillating flow about 
circular cylinders. In addition to modifying the long standing Morison Equation, this 
data also provided extreme insight into transverse forces on the cylinder. Performed at 
Reynolds numbers of prototype magnitude and across a wide spectrum of Keulegan- 
Carpenter [Ref. 6] numbers, these experiments provided a data base to which all 
subsequent publications and literature on the subject has been based. Of particular 
significance and not predicted by the Morison Equation was the documentation of 
cycle-to-cycle variations in pressures and forces across the cylinder in pure harmonic 
flow (see Figure E.1). These variations, particularly significant in the transverse 
direction, were attributed primarily to limited coherence lengths and to a lesser extent 
wall reflections and boundary layers. 

The distorting effects of coherence length on the wake of a circular cylinder in 
uniform flow was clearly discussed by Gerlach [Ref. 7] and Graham [Ref. 8.] In 
oscillating flows, however, the von Karman street forms perpendicular or transverse to 
the flow direction in the range 8<K,<13. This makes experimentation to identify the 
coherence length extremely difficult. In an effort to specifically examine the effects of 
spanwise coherence, O'Keefe (Ref. 9] in 1986 undertook similar experiments to 
Sarpkava [Ref. 5] but with the minimum aspect ratio (L/D = 2.0) that could be 
accurately instrumented. These nearly two dimensional experiments provided 
extremely accurate data that was not blurred by the influence of a limited coherence 
length. Cycle-to-cycle variations in in-line and transverse forces remained though; 
particularly in the Keulegan-Carpenter number range 8< K, < 13. 

Knowing that at any given point in time the forces acting on the cylinder is 
governed primarily by one dominant vortex, a detailed investigation of one potential 
vortex in the vicinity of a circular cylinder in harmonically oscillating flow was dictated. 
After an investigation of the single vortex case the follow-up was to investigate two 


potential vortices in the vicinity of a circular cylinder. 
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B. POTENTIAL FLOW MODELING 
1. One Vortex Model 


In the complex plane the potential flow field for a vortex in the vicinity of a 


circular cylinder in harmonically oscillating flow is given by the following expression: 


? 


f(z) = AS EA@E NZ <) cs iscin(z apr AN iso in(z —C-) & 


Zo 


where: c - radius of the cylinder 
I’ - strength of the vortex 
@ - frequency of the oscillating flow 


Unax > Maximum amplitude of the flow velocity 
Z, - location of the vortex 


In nondimensional form this expression becomes: 


PS) = sin(2ney(G + 7) + iKin(g - C)-iKin(g- eC 
Co 
where: ie Weel.) 
Beat) 1 (T - period of the flow) 
C = 2c 


This expression can be differentiated with respect to ¢ to determine the flow 


field velocities. At = , the expression for the velocity of the vortex becomes: 





ce = u-iv = sin(Q2nty1--4,)- Ah (3.3) 
C= Co > ry 
where r, is the radius between the real and image vortex and Is given by: 
V 0 an 
Co 
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(Note that the second term of eqn 3.2 is dropped since a vortex does not 
impart velocity to itself.) 


Calculations for the forces acting on the cylinder were made by integrating the 


Blasius Force Equation: 
7 = 2) edie -. O fF as 
D-ilL=— <-)-dz + ip—~¢9f dz 3:9 
IL = > pHa ae + ipa} - 


Nondimensionalized and evaluated at C=¢., this equation reduces to the 
following form: 








Dib’ = 2% cos(2nt) + BHR @ (¢ -1L) | (3.6) 
K. T : 


where: K. = (I We yVCZ¢) 


T 
“max 


2. Two Vortex Model 
Invoking the same procedure as the one vortex model it can be shown that the 


velocity of two vortices A and B in the vicinity of a the circular cylinder is given by: 


of = univ = sin(2nt)(1 — by — (1K _.) - (#8 _) + (8 ~) (3.7) 
oe C=O OA CAT Tan Ca Sp GAT High 
of =u-iv = sin(2t)(1 -—4 (SS Ee) (3.8) 
% C=C SRE, SaoSa SBE 


3. Rankine Vortex Model 
A standard rankine vortex model was used with linear velocity distribution 
between 0 and rog,,. The value of the core radius Was chosed to be 0.10 times the 


radius of the cylinder generating the vortices. 


4. Gaussian Time Dependent Vortex Model 
A standard Gaussian time dependent vortex model was used where the 


generalized expression for the velocity distribution is: 


r, 2 
= — ao aa a 
Vv = 2 it = exp( 5 \} (3.9) 





Here t. is not the clock time of the flow but the life time of the vortex. It has 
been shown that differentiating this with respect to the radius vields r= DR) (vt_) as 
the expression for the core radius of the vortex. Nondimensionalizing eqn 3.8 in terms 
of the oscillating flow parameters, and examing the velocity of a vortex from eqn 3.3. 


the velocity of a Gaussian time dependent vortex is given by: 


a. iv = sin(2ne)(1-—)~ FR (1 - exp( ~ MER) G10 
— = uiv = sin —l)- - exp “= 
aC _= Co Ly [oq tT) 





For the modeling the Reynolds Number was chosen to be Re= 10,000 and t, 
was set such that the core of the vortex was r =0.10 at inception. 
5. Numerical Integration and Accuracy 
In this analysis a finite differencing technique was employed. The study of 
potentially chaotic systems requires a high degree of accuracy. Due to its guaranteed 
stability and overall efficiency, the Adams-Moulton predictor corrector method was 


selected to convect the vortices: 





Znp = Zn 7 (99 In 7 99 Gn-2 F 374.3 7 7 Gna) 7 (3.11) 
z=2,,+(9q,. + 19 — 5 eacieyo oe (3.12) 
n n-1 Unp Onel Gn-2 * 4n-3 54 ale 


In addition to ensured stability, this method provided fourth order accuracy 
with local error O(h>). With 1440 time steps per cycle, the local error remained at the 
limits of double precission machine accuracy. Initial startup was made using the fourth 
order Runge-Kutta scheme. All continuations were based on the last four calculated 


locations of the vortices, therebv ensuring accuracy and duplication at any cycle. 


Vy 


While discussing accuracy it must be noted, as mentioned earlier, that by 
definition chaotic behavior 1s non-integrable. Whatever scheme used can only illustrate 
qualitatively the onset of chaos. If used for comparison universally, all methods yield 
equally valid qualitative results. The effort devoted here to accuracy is therefore to 


ensure the validity of the results and add to overall integrity of this report. 
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IV. RESULTS 


eee GENERAL 

From Sarpkaya [Ref. 5] and other studies, it is concluded that the region of 
special interest in the field of oscillating flows is confined to $< K.< 13, an interval in 
which a half transverse vortex street is generated. It is because of this reason that in 
the present numerical experimentation K, was chosen equal to 10.0. 

The vortex path is a function of both the imposed ambient flow field and the 
velocities induced by the vortices (real or image). This coupling was first investigated 
through the use of a single vortex in uniform flow in the vicinity of a circular cylinder. 
Figure E.2 illustrates the path of the vortex for various vortex strengths (K = 0.5, 0.625, 
0.75). It is seen that the capture oy the image vortex (partial looping) lasts onlv until 
the flow field around the cvlinder pushes the vortex away from the image. It is also 
seen that this capture is a function of the initial location of the vortex. When the 
vortex is captured, sharp cusps appear where the direction of motion of the vortex is 
quickly reversed. These characteristics played even more decisive roles in the unsteady 
case. With the exception of the foregoing, the vortex strengths used in the numerical 


experimentation were all set at K=0.50. 


B. ONE VORTEX MODEL 

Initial experimentation included varying the starting locations and times of the 
single point vortex in order to observe the vortex path and force coefficients. Figures 
E.3 through E.6 show the force coefficients for various starting locations and times and 
varying numbers of cycles. These figures clearly show cycle-to-cycle variations in both 
the inline and transverse force coefficients. In comparison to experimental force traces 
shown in Figure E.1, of comparable K., the variations are of similar magnitude and 
shape. 

Understanding the origin of these cycle-to-cycle variations required a scrutiny of 
the vortex paths themselves. Figures E.7 and E.8 show two starting positions of the 
vortex of (0, 1.5083) and (0, 1.5084). The oscillating flow velocity was set at U, = or 
(superscript indicates the initial x-direction). Figures E.7(a) and E.8(a) show that for at 
least the first five cycles of the flow the vortex paths are virtually identical. Also noted 
in this figure are sharp cusps where the overall direction of the vortices is abruptly 


changed. 


Zl 


Although the difference in the begining locations of the sixth cycle is virtually 
inperceptible {C=(0.995, 0.723) vice ¢=(0.970, 0.753)}, the vortex in Figure E.7(b) 
forms a sharp cusp at the top of the cylinder (geographical location of the maximum 
velocity at any time in the flow) but the vortex in Figure E.8(b) is captured by its 
image and orbits the cylinder. The ending points for cycle six are grossly different {see 
Figures E.7(c) and E.8(c)} and the vortex paths never again become similar. The finite 
perturbation in starting location has vielded explicitly different results. Nevertheless, 
the vortex paths in Figures E.7(c) and E.8(c) remain bounded. Extending the 
calculations to 100 cycles of the flow, Figure E.9 shows that while neither vortex paths 
are identical, both remain bounded. They clearly seem to slide in and out of chaotic 
behavior, preceded and followed by periods of stability. 

Figure E.10 is a plot of the magnitude of the velocity of the vortices. It is again 
evident that these two vortices have virtually identical velocity magnitudes until the 
sixth cycle of flow. The capture and looping depicted in Figures E.7(b) and E.8(b) 1s 
also evident in Figure E.10. Breaking down this velocity magnitude into its u and v 
components, Figures E.11 through E.13 show that overall periods of stability are 
intermixed with these vortex captures. 

Figures E.14 through E.19 and E.20 through E.25 are similar investigations that 
document two starting points of the vortices where both initial x-locations were set at 
ands Sa 


X max’ 
respectively. Figures E.26 through E.31 are similar investigations that document two 


x=0.0 and the oscillating flow velocity was set at U..= Uae 


Starting points of the vortex in the vicinity of x= 1.0, y=1.0, and the oscillating flow 
velocity was set at U,=07. These investigations support the conclusions of the first 
investigation and demonstrate once again a finite precission input does not yield a 
finite precission output. They also demonstrate the periods of stability that precede 
and follow chaotic behavior, (i.e., chaos in order and order in chaos). 

As seen in the proceeding set of experiments, the looping or continued looping of 
the vortex tightly around the cylinder, was a function not only of the location of the 
vortex but also of the current state of the oscillating flow. This suggests that at a 
given initial location, a variance of initial velocity or phase shift of the oscillating flow 
could drive the system from stable to chaotic behavior at a given point in time. Figure 
E.32 is a series of vortex paths that supports this concept. Four vortex paths were 
generated with the initial location of x=y=1.0 and a variance of initial velocity by a 


phase shift fromu@™ S50 ae atemiev (U,=07). It is clear that while the Figure 


a: 


E.32(a) documents stability even after ten cycles, the cusps and asymmetric tracks in 
Figure E.32(d) demonstrate the mixture of stable and chaotic regions. 

Figure E.33 depicts the results of a similar investigation where the initial location 
of the vortices is at (0, 2.2192) and (0, 2.2193). These results also show that the 
chaotic behavior depends not only on the vortex location but also on the flow 
environment. In fact, the two are coupled inseparably. 

The stability and svmmetrv of Figures E.32(a) and E.33(a) also suggest, however, 
that chaotic motion may not be inevitable. They infer that the location of the vortex 
and field flow may be tuned in such a way as to permanently avoid chaos. This would 
be paramount to permanently balancing on a razor’s edge. However, such a balance is 
unstable. As seen in the first investigation, the slightest disturbance is sufficient to 
push the system into chaotic behavior. 

The final phase of the one vortex experiment was to examine the effects of higher 
order models of a potential vortex in an effort to better understand vortices in nature. 
Figure E.34 is a comparison of the results obtained with ideal, Rankine, and Gaussian 
time dependent vortex models. It is clear from these figures that while the onset of 
chaos may be different, the chaotic behavior and susceptibility to chaos remains 


virtually the same. 


C. TWO VORTEX MODEL 

Separated flow past a bluff body generates pairs of counterrotating vortices. In 
both steady and oscillating flow about a circular cylinder, the first vortex pair is formed 
symmetrically and is significantly stronger than subsequent asymmetric vortices. 
Additionally, it has been documented that in oscillating flow in the region 8<K,< 13, 
a half von-Karman vortex street forms perpendicular or transverse to the flow 
direction. In this section of the experimentation, two vortices were placed Woh Weve 
Vicinity of a circular cylinder immersed in a harmonically oscillating flow. The vortices 
Were tracked through various numbers of cycles of flow. The locations and strengths 
of the vortices were varied in order to examine the mechanism for the formation of a 
half von-Karman street. Additionally, the stability of the vortices as discussed in the 
one vortex model was examined. For all of the investigations, the characteristic of the 
flow was set at K,=10.0. The strength of the vortices was set at K=0.50+6 
(although contrary in direction), and their locations at Ga = (a.a) and 
Cp =(at6,-at6), where 0.80<a<1.414, and 6< <1. Based on flow visualization, the 


untiall velocity was set at U.= U__.,. 
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The first part of the investigation required the analysis of vortices, symmetric in 
strength and location. Figure E.35 demonstrates that within the regions examined, the 
mutual induction velocities between the real vortices dominate the advection of the 
vortices. It is also clear that the vortices are quickly convected away in parallel paths. 
The initial radial location determunes the separation distance of the two parallel paths. 
This separation distance conibined with the strength of the vortices determines the rate 
at which the vortices are swept away. 

Figures E.36 through E.38 are of a series of investigations where the vortices 
have equal strengths (K,=Kp) but the locations are asymmetric (|§pl=€, X |Sa]). 
The overall result is a general curving of the vortex path toward the vortex that had 
the larger initial radial distance from the cylinder. This is a result of the image vortices 
having a greater effect on the vortex that was closer to the cylinder. It is seen from 
Figure E.37 that a 10% difference (€ =0.90) in the two starting locations is sufficient 
to significantly curve the vortex paths. The vortex paths do not gravitate to each other 
though and the vortices do not become coincident. 

Figures E.39 through E.42 are of a series of investigations where the vortices 
have symmetric initial locations ([,4 =p) and asymmetric strengths (Kp=ct, x Ky). 
The result here is a curvature similar to that of the previous investigation. From . 
Figure E.41 it is seen that only a 5% difference (€,=0.95) in the vortex strengths is 
sufficient to significantly alter the vortex paths. Clearly the results indicate that the 
system is much more sensitive to asymmetry in vortex strengths than to asymmetry in 
vortex location. Again the vortices do not become coincident. 

Figure E.43 is a combination of the previous investigations where there is 
asymmetry in both the vortex strengths (Kp=&.*Ka) and in initial location 
(ISpl= EX Sal): Again both vortex tracks bend toward the dominant vortex but as 
expected, with a more severe curvature. Although in real flow two vortices are 
generated every time the oscillating flow direction changes, this study of two vortices 
moving together to either the top or bottom of the cylinder is particularly significant. 
It suggests the evolution of a half von-Karman vortex street transverse to the flow 
based on the asymmetry of the vortices in the street. 

Finally, an investigation was made where the asymmetry of the vortices was so 
extensive that the vortex paths curved back to the vicinity of the cylinder. Although in 
viscous flow the life of a vortex is finite, combinations of additional vortices coupled 


with the flow field make this situation plausible and warrants its investigation. Figure 
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E.44 documents the sharp cusps, capture, and chaotic behavior of the vortices as 
discussed in the previous section. With the additional vortex (and image) the coupling 
is much more complicated between the geographical location of the vortex and the 


flow field. The identification of chaotic behavior is none the less possible. 


aS 


(2 


V¥. CONCLUSIONS 


The results presented herein warranted the following conclusions: 


Finite precision information about the initial characteristics of a vortex (or of 
larger number of vortices) provides no finite-precision information about the 
characteristics of flow at later times and that the lack of spanwise coherence of 
the vortices is not the sole cause of cycle-to-cycle variations in the measured 
quantities. | 

An idealized vortex model suggests that the sensitivity of the characteristics of 
the flow to random disturbances superimposed on the motion of vortices at an 
earlier time is the main flow feature responsible for the observed cycle-to-cycle 
varlations. 

The chaotic behavior of the vortex in oscillating flow occurs in regions of 
susceptibility similar to those found by Dowell and Pezeshki [Ref. 4] with 
similar sensitivity. 

The susceptible regions for chaos in oscillating flow are coupled inseperably 
with the location of the vortices and the state of the environmental flow field. 
In other worlds, for a given flow phase there is a vortex location that will lead 
to chaos. Likewise, for every vortex position near the cylinder, there is a 
particular flow phase that will eventually lead to chaos. 

The formation of a half von-Karman vortex street in the range of 8<K_<13 
may be attributed to the asymmetry in the formation of the vortices. 

The extent to which the sensitivity of the flow to the conditions prevailing in 
the previous cycles is damped by finite core effects and viscosity has been 
investigated and it has been shown that they change only the time of onset of 


the chaos, not its invevitable occurrence. 
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APPENDIX A 
ONE VORTEX COMPUTER MODEL 


RAKRKARRAKKKKKAKKAKRKKKRKKKKAKKRKKKKKKRKKKAKKAKKKRKKKKKKAKKKKKKKKRKKKKKKKKRKKK 


* 


* WO Ret Eex 
* 


* BYoLT. Wl. MeEcoy 
* 


KAKKKKKKKKAKKRKKKKKAKKKKKKKKKKKKKKKKRKKKKKRKKAKRKKKKKKKARKKKKAKKKKKKAKAKKR 


DESCRIPTION 


THIS PROGRAM TRACKS A VORTEX IN COMPLEX SPACE AROUND A UNIT CYLINDER 


* 


* 


* 


* 


* 


* 


* 


IN OSCILLATING FLOW. INITIAL LOCATION, STARTING TIME, AND THE NUMBER 


OF CYCLES OF THE OSCILLATING FLOW ARE REQUIRED INPUT. 


PHaeepeeGkhia, 15 DESIGNED TO RUN ON DISSPLA BY USING THE FOLLOWING 
Pies be ease COMMAND: (DISSPLA VORTEX'. IT 1S CURRENTLY 
DrveNvstONeD TO RUN A MAX OF 10 CYCLES AT 1440 STEPS PER CYCLE. 
2.0M OF VIRTUAL MEMORY IS REQUIRED TO RUN THE PROGRAM AT MAXIMUM 
CAPACITY. THE PROGRAM IS INTERACTIVE FOR INPUT OPTIONS. 


PeieeoOtouailo ARE SET IN THE MAIN PROGRAM. PREVIOUS INPUT OPTIONS 
Peter ROM FILE 8FILE FIZ5>FOOl*.. FAILURE TO GENERATE THIS FILE 
PRIOR TO INITIAL RUNNING RESULTS IN A ERROR MESSAGE AND TERMINATION 
Cte PRUOR TOeRUNNING, PIEE “FILE FIZ5FO01" SHOULD BE 


INITIATED CONTAINING 4 COMPLEX NUMBERS, 1 REAL NUMBER, AND 1 INTEGER. 


Crete es INPUL THE PROGRAM DOES NOT NESD TO BE TO BE RECOMPILED 
FOR EACH RUN. 


TolomepechaAM ITS WRITTEN IN DOUBLE PRECISION AND EMPLOYS THE MILNE 
People ton-CORRECTOR CONVECTION SCHEME, WITH INITIAL POINTS GENER- 
ATED BY A 4TH ORDER RUNGE-KUTTA SCHEME. 


Pail 


Mm © OF OY ©). Ty OF OF OO O° 20 Bere 


an Gs-0-O © 


Om 
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OUTPUT OF THE PROGRAM INCLUDES THE LAST FOUR VALUES OF THE LOCATION 
OF THE VORTEX, FINISH TIME, AND THE NUMBER OF CYCLES. [HE sEeVAEua> 
ARE WRITTEN IN FILE 'FILE FT25FO01' AND ALLOW FOR OBSERVATIONS BEYOND 
THE MAX LIMIT OF 10 CYCLES BY CONTINUING WITH ADDITIONAL RUNS. 


OUTPUT FOR THE VORTEX PATH AS WELL AS FORCE COEFFICIENTS AND 
VELOCITY MUNIPULATIONS ARE MADE GRAPHICALLY. SELECTIONS OF PLOTS 


FOR OUTPUT ARE MADE INTERACTIVELY. PLOTS MAY BE MADE ON 
A TEK618 OR SENT TO A METAFILE (PRE-EDIT SUBROUTINE 'GRAPH'). 


ALPHABETICAL LISTING OF VARIABLES 


VARIABLE (SUBROUTINE USED) DESCRIPTION. 


oe 028 8@ @ &@ @ @ 8 @ een ee 82 @@ @ @ HO O@OmUcUGHhUM OCU HOC HmLhLC HO FT FC HF HF FHF FFF HF FF HF HF GF HF He HH HeULDODmDUDLUOUDUCMDLCCHDLCOCDCUC HhmUMHMhUCUM HhUCUMAMrhCUcOrhCUCUchFOrhUCUchMhCUc hUCUchFOhUchFOhUchhUCUchUchhUchhUh}hUh!F 


CDG. ctetnemenene (GRAPH)NON-DIMENSIONAL ARRAY OF INLINE (DRAG) 
FORCE COEFFICIENTS FOR VORTEX POSITION Z(N). 


CIOGND sagem (GRAPH)NON-DIMENSIONAL ARRAY OF PERPENDICTULAR 
(LIFT) FORCE COEFFICIENTS FOR VORTEX POSITION 
Z(N). 

CORE sn is. sietete ss (FORCE ,GRAPH)COMPLEX VARIABLE CONTAINING AS ITS 


PARTS THE NON-DIMENSTONAL FORGE: GCErriGiaiilos. 


DELTAS. o.ane (STEP ,Q,OUTPUT, FORCE ,GRAPH) INCREMENTAL STEP 
IN NON-DIMENSIONAL TIME (INVERSE OF STEPS PER 
CYCLE). 

DELZ.........(STEP)COMPLEX INCREMENTAL STEP OF THE LOCATION 
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OF THE VORTEX {PROM Zi=TO 22). 


DFDZ........-(Q)COMPLEX PARTIAL DERIVATIVE OF THE COMPLEX 
PUNCRIONSE WITH TO RESPECT Z. 


DESPA Nave cones (STEP ,GRAPH)ARRAY CONTAINING THE TOTAL DISPLACEMENT 
OF THE VORTEX FOR A CORRESPONDING LOCATION Z(N). 


erect 6s ss (Q,FORCE)CONSTANT COMPLEX NUMBER = (0.0,1.0) 


JESS ae (a eens (GRAPH)LEFT HAND VALUE OF THE TIME AXIS OF THE 
PORGE GCOERPICIEND PLOT. 


JUSS) Sg Geneeponeerae (GRAPH)RIGHT HAND VALUE OF THE TIME AXIS OF THE 
FORCE CORRS IGIENT PLOT. 


ISTEPS.......(ALL)NUMBER Oe eRiGRrUoiinleolero PER CYCLE TO BE 
TAKEN FOR ALL ITERATIONS. 


Kie4... ee rar (RUNGE) RUNGE-KUTTA CONSTANTS 

RAGA. oss. .(ALL)STRENGTH OF THE VORTEX (GAMMA/(2*PI*UMAX) ) 

OMe ss cash (ALL) KEULLEGAN-CARPENTER NUMBER FOR THE OSCILLATING 
FLOW. 

i ee (INPUT,GRAPH) LOGICAL OPERATOR. 

NUMGYE....... (INPUT) TOTAL NUMBER OF CYCLES TO BE EVALUATED. 

NUMPTS....... (ALL) TOTAL NUMBER OF POINTS TO BE DETERMINED. 


(NUMCYC * ISTEPS) 


Promotes. + so (GRPAH)INITIAL PHASE SHIFT OF THE FLOW (0<PHSHFT<1) 


Gi) 7GY “OCF - Gy “Gy? (Ey. eGle een Ol Cy. Co et Cy Ot) Om) Ole CY em 20. SO): Oe CC COA 8) ~~ OD) SOI Oe SO) Os ON Sea) eee 


Pils ee ee (ALL)CONSTANT = 3.141592654 


0. eaten ae (Q)FUNCTION SUBROUTINE WHERE Q IS THE COMPLEX 
VELOCITY OF THE VORTEX. 


OMBGCN jis ncuane nce (GRAPH) ARRAY CONTAINING THE MAGNITUDE OF THE 
COMPLEAS VELOCITY OV cena 4 


OV EIGIION,) .acucues (STEP ,GRAPH)COMPLEX ARRAY CONTAINING THE VELOCITY 
OF THE VORTEX AT ANY POINT Z(N). 


SWirere eG banc (Q) VECTOR BETWEEN THE VORTEX AND ITS IMAGE 
WITHIN TRE CY EINER 


UP oo Oe (GRAPH)NON-DIMENSIONAL TIME (0-1 FOR A CYCLE). 
TE. eee eee eee (RUNGE ,Q)NON-DIMENSIONAL ENVIRONMENTAL FLOW TIME. 
bse ue peimen eters (STEP ,Q,FORCE)NON-DIMENSIONAL TIME CHANGED TO 


RADIANS FOR COS/SIN EVALUATIONS (0-2*PI FOR A CYCLE). 


) TT ie pie re cae (OUTPUT )NON-DIMENSIONAL FINISH TIME OF THE RUNS. 


StS coker © .. (ALL) INITIAL NON-DIMENSIONAL STARTING TIME FOR THE 
OSCILLATING FLOW. 


IS 6 GS Coc (Q)NON-DIMENSIONAL TIME LIFE OF THE VORTEX. 
(TE + CONSTANT FOR R-CORE AT INCEPTION) 


Oe ecopere cries vers (GRAPH) X-DIRECTION VELOCITY MAGNITUDE. 
Vice cesses. « (GRAPH) Y-DIRECTION VSLOGCITY ei Gianuber 
Divs atiaiea cone nate ... (ALL) INSTANTANEOUS COMPLEX LOCATION OF THE VORTEX. 
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OD 


Nya Aa Aa YY 


O 


ZN ane a ere (GRAPH) IMAGINARY(Y) LOCATIONOF THE VORTEX. 


AS. (Q)COMPLEX PREDICTED VALUE FOR Z(N) 
ZU, on cst (GRAPH)REAL(X) LOCATION OF THE VORTEX. 
Fyre Ree ass (Q)DUMMY ARGUMENT FOR THE FUNCTION SUBROUTINE Q 


WHICH CONTAINS THE PASSED VALUE OF Z(N). 
(RUNGE) DUMMY ARGUMENT CALL FUNCTION SUBROUITNE Q 


a MAIN PROGRAM a 
KREAEKKKAKKKKKKKKKKKAKKAKKAKKAAKKAKKKKAKKKAKAKKKKAKKKAAKKAKKRAAKAAKKKKRAKKRKK 


IMPLICIT REAL*8 (A-H,0-Z) 

REAL*8 KC,KAPPA 

REAL*4 DISP(14409) ,CD(14409) ,CL(14409) 

COMPLEX*16 Z(14409) ,QVECT(14409) . 

COMMON /FIRST/ TTI,ISTEPS /SECOND/ NUMPTS /THIRD/ KC,KAPPA 
l /FOURTH/ PI 


SEe CONSTANTS 


KAPPA = 0.5D0 

Reva o. DO 

PI = 3.141592654D0 
ISTEPS = 1440 


10 CALL INPUT (Z) 
CALL STEP (Z,QVECT,DISP) 
CALL OUTPUT(Z) 
CALL FORCE(Z,CD,CL) 
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CALL GRAPH (Z,QVECT,DISP,CD,CL) 


WRITE (6,*)'DO YOU WISH TO RE-RUN THE PROGRAM? 1-YES, O-NO! 
READ (5,%*) Ll 
IF (L1.NE JR eeonorzo 
WRITE (6,*) ‘EXECUTION CONTINUES ...... 
GOTO 10 
20 CONTINUE 


CALL DONEPL 
Sor 
END 


KKEKKKARKKKKKKKAKKAKKAKRKAKARKAKRKKRKRAKRAKRKKKRRRAKRKRKRKKRRKRKKRKKKKKKKKKKKKKKKAKKRKKAK 


* SUBROUTINE INPUT * 


KAARKKRKKKKKKKRKKKKKKRKAKRKAKKARAKKKRKAKRARAKKAKKAKAKKKKAKAKRKKKKRKAKKKRKKKKKKKAKRKAKE 


THIS SUBROUTINE INTERACTIVELY INPUIS THE VALUES FORM@PHE Tires 
LOCATION OF THE VORTEX, STARTING TIME, AND THE NUMBER OF CYCLES 
TO BE EXAMINED. 


SUBROUTINE INPUT (Z) 

IMPLICIT REAL*8 (A-H,0-Z) 

COMPLEX*16 Z(14409) 

COMMON /FIRST/ TTI,ISTEPS /SECOND/ NUMPTS 


READ THE LAST 4 POINTS FOR CONTINUATION OPTION 


REWIND 25 
DO 10 N=1,4 

10 READ(25,*) Z(N) 
READ(25,*) TTI 
READ(25,*) NUMCYC 


a 


C CONTINUATION OPTION 


© 
WRITE(6,*) 'DO YOU WISH TO CONTINUE THIS RUN FORM THE LAST? 1-YES 
# Q-NO! 
READ(5,%*) Ll 
IF (L1.NE.0) GOTO 20 
WRITE(6,*)' INPUT THE NEW STARTING POINT Z(1): (REAL,IMAG)' 
READ(5,*) Z(1) 
WRITE(6,*)' INPUT THE NEW VALUE OF TTI:'! 
READ(5,*) TTI 
C 


C CALL SUBROUTINE RUNGE TO GENERATE THE FIRST 4 POINTS 
@ 
CALL RUNGE (2) 


C 
C INPUT THE NUMBER OF CYCLES TO BE ANALYZED 
20 WRITE(6,*)'THE LAST VALUE FOR NUMCYC =',NUMCYC 
WRITE(6,*)'DO YOU WISH TO CHANGE IT? 1-YES O-NO! 
READ(5,*) Ll 
IF(L1.NE.1) GOTO 30 
WRITE(6,*) ‘INPUT THE NEW VALUE OF NUMCYC:'! 
READ(5,*) NUMCYC 
30 CONTINUE 
C 
NUMPTS = NUMCYC * ISTEPS + 4 
C 
RETURN 
END 


KRRAKKKKRKRKRKKRKRKARKRKKKRRRKRRRRKRKRRKRKRKRKKRKRRRKRRKRRRKKRRKKRRRKAKRKRKRRAKRKKRKKRRRKKAKRKKA 


* SUBROUTINE RUNGE S 


KRAKKAKKKKKKAKKAKKKAKEKKAAKRKAKKKAAKKRKRKAKKKKKKKARKAKAKRKAKKAKKRRKRKKRRKAKKKAKAKKRARKAARK 


Qi, (5 “Cy 10) 


C THIS SUBROUTINE USES A 4TH ORDER RUNGE-KUTTA MEU OD 1ONGr ieee laine 
C FIRST FOUR VALUES7One4e 


€ 
SUBROUTINE RUNGE (Z) 
IMPLICIT REAL*8 (A-H,0-Z) 
COMPLEX*16 Z(14409),Q,K1,K2,K3,K4,ZP1,2ZZ 
COMMON /FIRST/ TTI,ISTEPS 
€ 
DELTAT = 1.D0/DFLOAT(ISTEPS) 
C 
C RUNGE-KUTTA EVALUATION FOR FIRST FOUR POINTS 
DO 5 N=1,3 
Ad AON 
TE] = (DELOAT(N=1) * bELaan oon) 
Kl = DLETAT * (Zz ememn) 
ZZ = Z(N)+0.5D0*K1 
TE2 = TE1+ 0.5DO*DELTAT 
K2 = DELTAT 4uG(277Ee 
ZZ = Z(N)+0.5D0*K2 
KS = DELTAT “90122 2 Ge2) 
ZZ = Z(N) + K3 
TE4 = TEl + DELTAT 
K4 = DELEAT oO ez ere) 
C 3 
5 Z(N+1) = Z(N) + (K1-2.D0*K2+2.DO*K3+K4)/6.D0 
C 
C 
C 
RETURN 
END 
C 
€ 


C RRREAKKKKKRKAKKAKKKKAKAKRKAKRKKKAKKRKKKAARKKKRKAKRKRKKKAKKKRARAKKKAKKRRKRKKKKAKK 


Grex SUBROUTINE STEP * 
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THIS SUBROUTINE PROGRESSIVELY CALCULATES THE VALUES OF Z(N), BASED ON 
PeoReenerOR=CORRECTOR CONVECTION SCHEME. 

iY BESO DETERMINeSS THE VELOCITY VECTOR, MAGNITUDE AND PARTS FOR 
FUTURE MUNIPULATION AS WELL AS OVERALL DISPLACEMENT MAGNITUDES. 


SUBROUTINE STEP (Z,QVECT,DISP) 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL*4 DISP(14409) 

REAL*8 KC,KAPPA 

COMPLEX*16 Z(14409) ,QVECT (14409) ,DELZ,ZP,Q,2ZZ 

COMMON /FIRST/ TTI,ISTEPS /SECOND/ NUMPTS/THIRD/ KC,KAPPA 


DELTAT = 1.D0/DFLOAT(ISTEPS) 
EVALUATE VELOCITIES AND DISPLACEMENTS FOR FIRST 4 POINTS 
pasri(1) = 0.0 
pomyen = 1,4 
TE = (DFLOAT(N-1) * DELTAT + TTI) 
ZZ Z(N) 


OVECT(N) = 0(2Z,TE) 
7 DISP(N+1) = 2(N+1) - 2(N) 


De 10 N = 5, NUMPTS 


iat = N-1 
NM2Z = N-2 
NM3 = N-3 
NM4 = N-4 


TE = (DFLOAT(N-1)*DELTAT+TTI) 


CALCULATE PREDICTED VALUE OF Z(N) 


oy) 


O 


Q 


ZP = Z(NM1)+(55.DO*OVECT GQIM1) =59- DO“OVEGT iZ) “ta7 900 orem 112) 
1 -9.DO*QVECT(NM4) )*DELTAT* 2.D0 * KC / 24.D0 


CALCULATE ACTUAL VALUE OF Z(N) 
DELZ = (9.D0*Q(ZP,TE) + 19.DO*QVECT(NM1) -5.DO*QVECT(NM2) 
1 +QVECT(NM3) )*DELTAT * 2.D0 * KC / 24.D0 
Z(N) = Z(NM1) + DELZ 
EVALUATE VELOCITIES AT Z(N) 
QVECT(N) = Q(2Z(N),TE) 


10 DISP(N) = DISP(NM1) + ABS(DELZ) 


RETURN 
END 


THIS FUNCTION SUBPORGRAM IS USED TO EVALUATE VALUES OF 0 


COMPLEX FUNCTION Q*16(2ZZ,TE) 

IMPLICIT REAL*8 (A-H,0O-Z) 

REAL*8 KC, KAPPA 

COMPLEX*16 ZZ, DFDZ, I,RV 

COMMON /FIRST/ TII,ISTEPS /THIRD/ KC,KAPPA /FOURTH/ PI 


ie = (cy jovel aL .ioval 


DELTAT = 1.D0/DFLOAT(ISTEPS} 
TT>°= Tene pexern 
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RV = 2Z - 1.D0 / DCONJG(2ZZ) 
DFDZ =DSIN(TT)*(1.D0-1.D0/ (ZZ*ZZ) ) -I*KAPPA/RV 


Q = DCONJG(DFDZ) 


RETURN 

END 
KREKKKKKKKKKKKKKKKKKRRKAKKKAKKAKKKKKKKKKKKKKKKKKRKKKRKKRKRKKRKKKRKRKRKRKRKKRKKKRKKKKAK 
7 SUBROUTINE OUTPUT x 
REKKAKKKKKKKKKKKKRKKKRKKKKKKRKKKKKKKKKKKKKKKRKKAKRKRKRKRKRKRKRKRARRKRKRRRKRKKRKKKKK 
THIS SUBROUTINE IS USED TO SAVE THE FINAL VALUES OF Z,TT, AND NUMCYC 
FOR POSSIBLE CONTINUATION OF THE PLOTS. 
THE VALUES ARE SAVED IN FILE 'FILE FT25FOO1'. 


SUBROUTINE OUTPUT (Z) 

IMPLICIT REAL*8 (A-H, O-Z) 

COMPLEX*16 Z(14409) 

GOMMen /FIRST/ TTI,ISTEPS /SECOND/ NUMPTS 


DabgaAr = 1,.D0/FLOAT(ISTEPS) 

NN = NUMPTS - 4 

TTF = DFLOAT(NN)* DELTAT + TTI 
NUMCYC = NN / ISTEPS 


REWIND 25 
DO 10 N=1,4 
NN = NN + 1 

10 WRITE(25,*) Z(NN) 
WRITE(25,*) TTF 
‘WRITE(25,*) NUMCYC 


SH 


7-7 AFA A AA A A TY 


RETURN 
END 


KREKKKAKKRKKRKRRKRARKRKRRKRRRKARAKRRRARRRRRRRRKRRRARKARRRRRKRRRRRRRRRRARRRRRRRRRRARKA 


x SUBROUTINE FORCE * 
KEKKKKRRRRKKRRRKKRKRRKRKRKREKKRKRKRKKRKRRRR RRR RRA RK KRKRKRKKRKRKKKRKRRRRRKRKRRRRERER 


THIS SUBROUTINE IS USED TO CALCULATE THE FORGE COEPRPICIENTS CD?) arp 
CL FOR THE CYLINDER BASED ON THE LOCATIONS OF THE VORTEX AND A 
FINITE DIFFEKRENCING "SCHEME. 


SUBROUTINE FORCE(Z,CD,CL) 

IMPLICIT REAL*8 (A-H, 0-Z) 

REAL*4 CD(14409) ,CL(14409) 

REAL*8 KC, KAPPA 

COMPLEX*8 COEF 

COMPLEX*16 Z(14409), I 

COMMON /FIRST/ TTI,ISTEPS /SECOND/ NUMPTS /THIRD/ KC,KAPPA 
1 /FOURTH/ PI 


Tl = (0 DO elo 
DELTAT = 1.D0/DFLOAT(ISTEPS) 


DO 10 N = 2, NUMPTS 


TT = (DFLOAT(N-1) * DELTAT + TTI) * 2ep0me=erd 


COEF = (2.DO*PI*PI/KC)*DCOS (TT) 
1 + ( (PI“KAPPA)/  Ke- Dera ee 
2 ((Z(N)-1.D0/DCONJG(Z(N)))- 
3 (Z(N-1) - 1.DO0/DCONJG(Z(N-1)))) 


COEF = CONJG(COEF) 
CD(N) = REAL(COEF) 
10 CL(N) = IMAG(COEF) 
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RETURN 
END 


RAAARAKAKKKAKKARARKKAKRAKRAKRARKRKKRKKARKAKAARAKRRAKAKAAKARAKKKKAKRARRAARKAKAKRERA 


x SUBROUTINE GRAPH . * 


Qo eh WE te AEE ee eT DOD eee, Te a te eee all el Ue Seg SO Oe, SS eee (OM eR TOR’ ee eR ae ean JO] Se Sen J (ee Se Oe 


SUBROUTINE GRAPH (Z,QVECT,DISP,CD,CL) 

IMPLICIT REAL*8 (A-H, 0-2) 

REAL*4 ZR(14409) ,ZI (14409) ,X(361) ,Y(361) ,QMAG(14409) 
REAL*4 DISP(14409) ,CD(14409) ,CL(14409) ,T(14409) 
REAL*4 U(14409) ,V(14409) , GARRAY(1) 

COMPLEX*16 Z(14409) ,QVECT(14409) 

COMMON /FIRST/ TTI,ISTEPS /SECOND/ NUMPTS /FOURTH/ PI 


DEFINE CYLINDER OF RADIUS 1.0 


DO 10 N = 1,361 
RAD = FLOAT(N-1) * PI /180.0 
X(N) = COS(RAD) 

10 Y(N) = SIN(RAD) 


DEFINE REAL SINGLE PRECISION NUMBERS FOR PLOT ROUTINES 


DELTAT = 1.D0/DFLOAT(ISTEPS) 


DO 20 N =1, NUMPTS 
ZR(N) =DREAL(Z(N) ) 
ZI(N) =DIMAG(Z(N)) 
QMAG(N) = ABS(QVECT(N)) 
U(N) = DREAL(QVECT(N)) 
20 V(N) = DIMAG(QVECT(N)) 


C 
c 


© iO: [Oly Ose Ore 
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GENERATE EEO@S 


CALE COMeRs 
CALL TEK618 


CALL BLOWUP(.3365655) 
CALL BLOWUP(.55) 
CALL BLOWUP(2.0) 


CALL NOBRDR 


VORTEX TRACK PLOT 


CALE GRAGeMCOMe) 


CALL NOCHEK 


CREE PAGE (1.0 ,8.5) 


CALL AREA2D(10., 5.71425) 


CALL YNAME('Y$',100) 
CALL XNAME('X$' ,100) 
CALL YAXANG(0.0) 


CALL INTAX&S 


CAEL KT ICKSCZ) 
CALL YII¢KsCZ) 


CALL GRAF(-7,1,7,-4,1,4) 


CALL CURVE(X,Y,361,0) 


GARRAY(1) = 


CALL SHADE(X,Y,361,45,GARRAY,1,0,0) 


CALL CURVE(ZR,ZI,NUMPTS,20000) 
CALL CURVE(ZR,ZI,NUMPTS ,0) 


CALL ENDPL(0) 


MAKE CONTINUATION SELECTIONS 


WRITE (6,%) 


CONTINUATION MEMO’ 
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WRITE(6, *) 


WRITE(6,*) 
WRITE(6,%*)' TERMINATION OF PROGRAM. ......cccceccecceces Q! 
WRITE(6,*)'RERUN PROGRAM......... NI ithe. oa 3: proses 1! 
WRITE(6,*)'CONTINUE WITH NORMAL SEQUENCE OF PLOTS..... zs 
WRITE(6,*) 


WRITE(6,*)'PLEASE ENTER YOUR SELECTION: (0,1,2)! 
READ(5,*)L1 

Peet. 1 )Goro 100 

iain EO.1)GOTOmIOO 

CONTINUE 


PHSHFT = TTI - DFLOAT(INT(TTI) ) 


DO 40 N = 1,NUMPTS 
40 T(N) = (DFLOAT(N-1) * DELTAT + TTI - PHSHFT) 


WRITE(6,*)'DO YOU WANT A FORCE PLOT? 1-YES, O-NO'! 
Reebok 
Deel NE “1)GOTO. 50 


P@nen, VS TIME PLOTS 


WRITE(6,*) ‘ENTER VALUE OF THE LEFT HAND SIDE OF THE TIME SCALE: ' 
READ(5,*) ILHS 

WRITE(6,*) 'ENTER VALUE OF THE RIGHT HAND SIDE OF THE TIME SCALE: ' 
READ(5,*) IRHS 


CALL AREA2D(10.0, 7.5) 
CALL YNAME('COEF$',100) 
CALL XNAME('T*$! ,100) 
CALL YAXANG(0.0) 

CALL INTAXS 


Gite) GRABGILHS,.5,IRHS,-4,1.,4.0) 
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CALL CURVE(T,CD,NUMPTS ,0) 
CALL DASH 

CALL CURVE(T,CL,NUMPTS ,0) 
CALL RESET ('DASH') 

CALL ENDPL(0) 


50 WRITE(6,*)'DO YOU WANT A Q-MAG VS DISP PLOT? 
READ(5,%*) Ll 
TF (Laie) GOne 60 


MAGNITUDE OF QO VS. DISPLACEMENT PLOT 


GhEL) ARBR20 (1070707 25) 
CALL YNAME('Qs' ,100) 
CALL XNAME('SS$',100) 
CALL YAXANG(0.0) 

CALL INTAXS 


CARLSGRMEr (2a) 
CHIL, GRAS OF 10 16G 4041s psn 


CALL CURVE(DISP,QMAG ,NUMPTS ,0) 
CALL ENDPL(0) 


60 WRITE(6,*)'DO YOU WANT A Q-MAG VS TIME PLOT? 
READ(5,*) Ll 
IF (L1.NE.1)GOTO 70 


MAGNITUDE Oi © Vs.) TIME Lor 


CALL AREA2D(10.0, 7.5) 
CALL YNAME('Q$',100) 
CALL XNAME('T*$' ,100) 
CALL YAXANG(0.0) 


1-YES, O=NGE 


ba OAS OL iC, 


2) 


CALL INTAXS 


Creer GRAM (O17 10,0,1.,20.0) 


CALL CURVE(T,QMAG,NUMPTS,0) 
CALL ENDPL(0) 


70 WRITE(6,*)'DO YOU WANT A U-VEL VS X-POS PLOT? 
READ(5,*) Li 
Pl oNe. 2 GOTO 80 


Wave ROCiTY POSITION PLOTS 


GALE AREAZD(10.0, 7.5) 

CALL YNAME('U-VELOCITY$! ,100) 
CALL XNAME('X-POSITIONS',100) 
CALL YAXANG(0.0) 

CALL INTAXS 


CAbEeGRAr(=7,1,7;710,2.,10.0) 


CALL CURVE(ZR,U,NUMPTS ,0) 
GAEL ENDPL(0) 


80 WRITE(6,*)'DO YOU WANT A V-VEL VS Y-POS PLOT? 
READ(5,*) Ll 
fe bl NS .1)GOTO 90 


VM=VBE@eCiiy POSITION PLOTS 


CALL AREA2D(10.0, 7.5) 

CALL YNAME('V-VELOCITY$' ,100) 
CALL KNAME('Y-POSITIONS! ,100) 
CALL YAXANG(0.0) 

CALL INTAXS 
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Toes, O-NO" 


Pexv85,, 


O-NO'! 


oe) 


CALL GRAR(=4,. 5307 =e) 


CALL CURVE(ZI,V,NUMPTS,0) 
CALL ENDPL(0) 


90 WRITE(6,*)'DO YOU WANT A U-VEL VS V-VEL PLOT? 1-YES, O-NO! 
READ(5,*) Li 
IF (L1.NE.1)GOTO 100 
U VS V VELOCITY PLOTS 
CALLIN ARSAZDN 7. 5507 15) 
CALL YNAME('V-VELOCITYS' ,100) 
CALL XNAME(!U-VELOCITYS$! ,100) 
CALL YAXANG(0.0) 
CALL INTAXS 


CALL GRAF(-10,2,10,-10,2.,10.0) 


CALL CURVE(U,V,NUMPTS ,0) 
CALL ENDPL(0) 


100 CONTINUE 


RETURN 
END 
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APPENDIX B 
RANKINE VORTEX COMPUTER MODEL 


Po oe UNC TION SUBROUTINE IS WRITTEN FOR USE WITH PROGRAM ‘VORTEX 
FORTRAN'. IT CONTAINS VELOCITY CALCULATIONS FOR A RANKINE 
WORTEX MODEL. IT MAY BE INSERTED IN PLACE OF THE IDEAL 

VORTEX MODEL (FUNCTION SUBROUITNE 9 IN 'VORTEX FORTRAN’ ) 


aoe ese FF eee SF SS FSS SF S22 2 SF SF SF ee ese Fees Few Fe VF ese FF FF FTFZ FZ SF FFT FT VZZFZwVZVwIZ_F2 2 S28 SF SF SF S28 S28 2 2 S28 S28 2 2 82 9 82 2 = = = 


a ae aseewreeww@t 2 2 2 2 2 2 2 2 2 2 2 2 es 2s See Sees 2 S2F es Fetes SF SF 2 SF SF SF SF SF Fetes 82 82 S222 FF FF SF ss SF Fw SF es Ses Ss Ss = = 


THIS FUNCTION SUBPORGRAM IS USED TO EVALUATE VALUES OF Q 


COMPLEX FUNCTION 0Q*16(ZZ,TE) 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL*8 KC, KAPPA 

COMPLEX*16 ZZ, DFDZ, I, RV, VINDUC 

COMMON /FIRST/ TTI,ISTEPS /THIRD/ KC,KAPPA /FOURTH/ PI 


RCORE = 0.1D0 
im—aOsD0),1.D0) 

DELTAT = 1.D0/DFLOAT(ISTEPS) 
foeeTe ~*~ 2,00 * PI 


Ea= 922) = D0 / BCONIG(ZZ) 
CHECK = ABS(RV) 


Pe eeReORe GT, CHECK) GOTO 10 


Vinibuc = =1,.D0 7 RV 
GOTO 20 


10 VINDUC = -RV / RCORE**2 


20 DFDZ =DSIN(TT)*(1.D0-1.D0/(ZZ*ZZ) )+KAPPA*IXVINDUC 


Q = DCONJG(DFDZ) 


RETURN 
END 
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APPENDIX C 
GAUSSIAN VORTEX COMPUTER MODEL 


iitos RUNCTION SUBROUTINE IS WRITTEN FOR USE WITH PROGRAM ‘VORTEX 

FORTRAN'. IT CONTAINS VELOCITY CALCULATIONS FOR A GAUSSIAN TIME 

PeeePeNTSVORTIEA MODEL. IT MAY BE INSERTED IN PLACE OF THE IDEAL 
VORTEX MODEL (FUNCTION SUBROUITNE Q IN 'VORTEX FORTRAN' ) 


THIS FUNCTION SUBPORGRAM IS USED TO EVALUATE VALUES OF Q 


COMPLEX FUNCTION 9*16(2Z,TE) 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL*8 KC, KAPPA 

COMPLEX*16 ZZ, DFDZ, I,RV 

COMMON /FIRST/ TTI,ISTEPS /THIRD/ KC, KAPPA /FOURTH/ PI 


e=CG1DG 1.0) 

DELTAT = 1.D0/DFLOAT(ISTEPS) 

Tr = TE * 2.D0* PI 

TVINV = 1.D0 /(TE + .497440685D0) 


Pue= 22 = 1.00 / DCONIG(ZZ) 
RVMAG = ABS(RV) 


ARG = -62.5DO*TVINV*RVMAG**2 
IF (ARG.LE.50.D0) GOTO 10 
ARG = 50.D0 


10 VMOD = 1.D0 - DEXP(ARG) 
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DFDZ =DSIN(TT)*(1.D0-1.D0/ (ZZ*ZZ) )-1* (KAPPA/RV) *VMOD 


Q = DCONJG(DFDZ) 


RETURN 
END 
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APPENDIX D 
TWO VORTEX COMPUTER MODEL 


KRRRKKARARAKRKKARKKKRRKKKKRRKRKARARKARAERKRKRKARRERERARRKERKERRARERKKRRKRRKRRRRKKEKE 
* 


* ; TW On VEO ok 
* 


a BY ei Wo  Meeoy, 
x 


EE ee SS te tet tte te ett ttt ttt tt tt. tet tt. tt.t8stitatatatetetotatatatetetetatetatetetatatatetatetetetetetas 


DESCRIPTION 


THIS PROGRAM TRACKS TWO VORTICIES IN COMPLEX SPACE AROUND A UNIT 
CYLINDER IN OSCILLATING FLOW. INITIAL LOCATION, STARTING TIME, AND 
THE NUMBER OF CYCLES OF THE OSCILLATING FLOW ARE REQUIRED INPUT. 


THE PROGRAM IS DESIGNED TO RUN ON DISSPLA BY USING THE FOLLOWING 
Diteseen nAkC COMMAND: 'DISSPLA TWOVOR’. IT IS CURRENTLY 
PEvehiStONvED TO RUN A MAX OF © CYCLES AT 1440 STEPS PER CYCLE. 
Z2.0M OF VIRTUAL MEMORY IS REQUIRED TO RUN THE PROGRAM AT MAXIMUM 
Ore eye) THE PROGRAM IS INTERACTIVE FOR INPUT OPTIONS. 


ALL CONSTANTS ARE SET IN THE MAIN PROGRAM. PREVIOUS INPUT OPTIONS 
Deena FILE ‘FILE FITS5R001'. FAILURE TO GENERATE THIS FILE 
PRIOR TO INITIAL RUNNING RESULTS IN A ERROR MESSAGE AND TERMINATION 
Cieeneee uN | PRIOR TO RUNNING, FILE “FILE FI35FO01" SHOULD BE 


* 


* 


* 


* 


* 


* 


* 


INITIATED CONTAINING 8 COMPLEX NUMBERS, 1 REAL NUMBER, AND 1 INTEGER. 


fete tahoe lNrUL THE PROGRAM DOES NOT NEED TO BE TO BE RECOMPILED 
FOR EACH RUN. 


THIS PROGRAM IS WRITTEN IN DOUBLE PRECISION AND EMPLOYS THE MILNE 
PREDICTOR-CORRECTOR CONVECTION SCHEME, WITH INITIAL POINTS GENER- 
ATED BY A 4TH ORDER RUNGE-KUTTA SCHEME. 
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‘@) 


a AO OQ MasG 7. +O) Oo “OF OA © © © (1) «Fy 


OUTPUT OF THE PROGRAM INCLUDES THE LAST FOUR VALUES OF THE LOCATION 
OF THE VORTICES, FINISH TIME, AND THE NUMBER @ OF CYCLES foto eer 
ARE WRITTEN IN FILE 'FILE FT35F001' AND ALLOW FOR OBSERVATIONS BEYOND 
THE MAX LIMIT OF 6 CYCLES BY CONTINUING WITH ADDITIONAL RUNS. 


OUTPUT FOR THE VORTEX PATHS ARE MADE GRAPHICALLY. PLOTS MAY BE MADE 
ON A TEK618 OR SENT TO A METAFILE (PRE-EDIT SUBROUTINE 'GRAPH' ). 


ALPHABETICAL LISTING OF VARIABLES 


VARIABLE (SUBROUTINE USED) DESCRIPTION. 

Dae een eters (STEP,Q,OUTPUT, FORCE ,GRAPH) INCREMENTAL STEP 
IN NON-DIMENSIONAL TIME (INVERSE OF STEPS PER 
CYGEE) 

DPD. ee (QA,QB)COMPLEX PARTIAL DERIVATIVE OF THE COMPLEX 


FUNCTION F WEIN TO Rest Sere. 


Nera eeo es can enero (QA,QB,FORCE)CONSTANT COMPLEX NUMBER = (0.0,1.0) 


PSA SRS sce (ALL)NUMBER OF INCREMENTAL STEPS PRR CYeCEEaro so 
TAKEN FORALL IITERATI Ons. 


Klass . 5 eee (RUNGE) RUNGE-KUTTA CONSTANTS 
KEEP eee see (ALL)STRENGTH OF THE VORTEX (GAMMA/(2*PI*UMAX) ) 
KAPPAB ae sneer (ALL)STRENGTH OF THE VORTEX (GAMMA/(2*PI*UMAX) ) 
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KC...........(ALL) KEULLEGAN-CARPENTER NUMBER FOR THE OSCILLATING 
FLOW. 

es crs Venemehens same (INPUT) LOGICAL OPERATOR. 

NUM CV Gis cis. a «16 (INPUT)TOTAL NUMBER OF CYCLES TO BE EVALUATED. 

NGMRAS . «ches (ALL)TOTAL NUMBER OF POINTS TO BE DETERMINED. 


(NUMCYC * ISTEPS) 
Oil, ee (ALL)CONSTANT = 3.141592654 


Otis sold ole (QA)FUNCTION SUBROUTINE WHERE Q IS THE COMPLEX 
VELOGITY OF THE VORTEX. 


ee raialietiew ssi’ (QB)FUNCTION SUBROUTINE WHERE Q IS THE COMPLEX 
VELOCITY OF THE VORTEX. 


OQVECTA(N)....(STEP,GRAPH)COMPLEX ARRAY CONTAINING THE VELOCITY 
OF THE VORTEX AT ANY POINT Z(N). 


OVECTA(N)....(STEP,GRAPH)COMPLEX ARRAY CONTAINING THE VELOCITY 
OF THE VORTEX AT ANY POINT Z(N). 


Ae Eevsa sie 6s 4.6 0-0 (RUNGE ,QA,QB)NON-DIMENSIONAL ENVIRONMENTAL FLOW TIME. 


Pitas) 56. 6. 0 5 9's (STEP ,QA,QB)NON-DIMENSIONAL TIME CHANGED TO 
RADIANS FOR COS/SIN EVALUATIONS (0-2*PI FOR A CYCLE). 


THR, - 0 ee ses (OUTPUT )NON-DIMENSIONAL FINISH TIME OF THE RUNS. 


0, (ALL)INITIAL NON-DIMENSIONAL STARTING TIME FOR THE 
OSCULEAT ING lon. 


é 7 Riad cae (ALL) INSTANTANEOUS COMPLEX LOCATION OF THE VORTEX. 
é . 

C ZB......+++..+(ALL) INSTANTANEOUS COMPLEX LOCATION OF THE VORTEX. 
C 

C ZIT coe (GRAPH) IMAGINARY(Y) LOCATIONOF THE VORTEX. 

C 

C A: EE (GRAPH) IMAGINARY (Y) LOCATIONOF THE VORTEX. 

€ 

e 2 (QA,QB)COMPLEX PREDICTED VALUE FOR Z(N) 

C 

C ZAR. . eee (GRAPH) REAL(X) LOCATION OF THE VORTEX. 

C 

C ZBRca ee eT ore (GRAPH)REAL(X) LOCATION OF THE VORTEX. 

© 

C Ze (QA,QB)DUMMY ARGUMENT FOR THE FUNCTION SUBROUTINE QA 
C AND QB WHICH CONTAINS THE PASSED VALUE OF Z(N). 

C (RUNGE) DUMMY ARGUMENT CALL FUNCTION SUBROUITNE QA 
CAND QB 

II DY Gee a ere pre ere) LtS 5 Goss o AOE 
é 

@ 

C 


CAKKARKKKKARKAKAKKKKAAKAKKRAAKKKKKAAKARKKKARARRKAKAKKAKKRRKKAKAKAKKAKRKAAKAAKRKERK 


ees MAIN PROGRAM as 
CARKARKKAAAKKAARARKAKKAKKARKRAAAKKKAKKKKAKRKAAAKKKKAKKARKARKKKKKAKKKKAKKAEKKKA RK 


é 
IMPLICIT REAL*8 (A-H,0-Z) 
REAL*8 KC,KAPPAA,KAPPAB 
COMPLEX*16 ZA(8649) , QVECTA(8649) 
COMPLEX*16 ZB(8649) , QVECTB(8649) 
COMMON /FIRST/ TTI,ISTEPS /SECOND/ NUMPTS /THIRD/ KC,KAPPAA, 
= KAPPAB/FOURTH/ PI 
C 
C SET CONSTANTS 
C 


Ci CY 02 OY OP 3” FE: oe 


KAPPAA = 0.5D0 

KC = 10.D0 

Hime 5. 14159265400 
towers = 1440 


10 CALL INPUT (ZA,ZB) 
CALL STEP (ZA,QVECTA, ZB,QVECTB) 
CALL OUTPUT(ZA,ZB) 
CALL GRAPH (ZA,QVECTA,2ZB,QVECTB) 


WRITE (6,*)'DO YOU WISH TO RE-RUN THE PROGRAM? 1-YES, 0-NO! 
READ (5,%) Ll 
IF (L1.NE.1) GOTO 20 
WRITE (6,*) 'EXECUTION CONTINUES ...... | 
GOTO 10 
20 CONTINUE 


CALL DONEPL 
SLOT 
END 


RRERKRKKAARRRARARKRRRKAKRRKRRKRKRKRKRRKARKRKRKRRRRRRRRKARKRKRRARRRKRKRKRRKRRRKRRRRKRKRRRRRRRKRRRRKREI 


~ SUBROUTINE INPUT * 


KRREARKKKKRKKKRKKKKARKKKRKKKRKAKKKRKRKRAKRAKRKAKRKKRRRKRRAKKRKERRARARKKRRKRKRRRRKRKRRRRRRRER 


THIiS@eeeeeUTINE INTERACTIVELY INPUTS THE VALUES FOR THE INITIAL 
Pee UGiMere THE VORTEX, STARTING TIME, ITERATION STEP SIZE, AND 
THE NWOMBER CF CYCLES TO BE EXAMINED. 


SUBROUTINE INPUT (ZA,ZB) 

MiP Eeer? REAL*8 (A-H,O-Z) 

REAL*8 KC, KAPPAA, KAPPAB 

COMPLEX*16 ZA(8649) ,ZB(8649) 

COMMON /FIRST/ TTI,ISTEPS /SECOND/ NUMPTS/THIRD/KC,KAPPAA, KAPPAB 
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READ LAST 4 POINTS OF EACH ARRAY FOR CONTINUATION OPTIONS 


REWIND 35 
Do 10 Nl Sie 

10 READ(35,*) ZA(N),ZB(N) 
READ(35,*) TTI 
READ(35,*) NUMCYC 


CONTINUATION OPTIONS 


WRITE(6,*) 'DO YOU WISH TO CONTINUE THIS RUN FORM THE LAST? 
+ Q-NO! 

READ(5,*) Ll 

IF (L1.NE.0) GOTO 20 


WRITE(6,*)'INPUT THE VALUE OF X : ZA(1)=(X,X)'! 
READ(5,*) X 

ZA(1) = DCMPLX(X,X) 

WRITE(6,*)'' INPUT THE VALUE OF EPSILON RADIUS: ! 
READ(5,*) ER 

X = ER * X 

ZB(1) = DCMPLX(X,-X) 

WRITE(6,*) ‘INPUT THE VALUE OF EPSILON KAPPA:' 
READ(5,*) EK 

KAPPAB = KAPPAA * EK 

S(O), Eiave 


CALL SUBROUTINE RUNGE TO GENERATE FIRST 4 POINIS POR BEACH ARRAY 


CALL RUNGE (ZA,ZB) 
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1-YEs 


C INPUT THE NUMBER OF CYCLES TO BE ANALYZED 
3 
20 WRITE(6,*)'THE LAST VALUE FOR NUMCYC =! ,NUMCYC 
WRITE(6,*)'DO YOU WISH TO CHANGE IT? 1-YES O-NO! 
READ(5,*) Li 
IF(LL.NE.1) GOTO 30 
WRITE(6,*) 'INPUT THE NEW VALUE OF NUMCYC:! 
READ(5,*) NUMCYC 
30 CONTINUE 


NUMPTS = NUMCYC ~*~ ISTEPS + 4 


RETURN 
END 


KRAKKAKAARKKKARKKKARKAKKAAKKARAKRARKAAARKRRARAARRKARRAKERKKARKKAARKKAKAKK 


* SUBROUTINE RUNGE | x 
KERR KKK RRR RRR RRR KERR RRR RRR RRR RRR KKK 


THIS SUBROUTINE USES A 4TH ORDER RUNGE-KUTTA METHOD TO GENERATE THE 
FPORSPSreuUR VALUES OF Z. 


CY “Os SG ey Cl ee Oy 0) 


SUBROUTINE RUNGE (ZA,ZB) 

IMPLICIT REAL*8 (A-H,0-Z) 

COMPLEX*16 ZA(8649) ,ZB(8649) ,QA,K1,K2,K3,K4,ZP1,22Z 
Corio FIRST/ TTL,ISTEPS 


C 
DELTAT = 1.D0/DFLOAT(ISTEPS) 
€ 
C RUNGE-KUTTA EVALUATION FOR FIRST FOUR POINTS 
c 


BOeseNa= 1,3 
“ZZ = ZA(N) 
ine = (PRLOnD (N-1) = DELDTAT 7.11) 
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Kl = DLETAT * QA(ZZ,ZB(N) , TEL) 
ZZ = ZA(N)+0.5D0*K1 
TE2 = TE1+ 0.5DO*DELTAT 


K2 = DELTAT * QA(ZZ,ZB(N) ,TE2) 
ZZ = ZA(N)+0.5D0*K2 

K3 = DELTAT * QA(ZZ,ZB(N), TE2) 
ZZ = ZA(N) + K3 


TE4 = TEL + DELTAT 
K4 = DELTAT * QA(2ZZ,ZB(N), TE4) 


ZA(N+1) = ZA(N) + (K1-2.DO*K2+2.D0*K3+K4) /6.D0 


YAY, =e RMON 

Kl = DLETAT * OB(ZA(N),2Z,TE1) 

ZZ = ZB(N)+0.5DO*K1 

K2 = DELTAT * OB(ZA(N),ZZ,TEZ) 

ZZ = ZB(N)+0.5D0*K2 

K3 = DELTAT * OB(ZA(N) (2zpenEz) 
ZZ = ZB(N) + K3 

TE4 = TEL + DELTAT 

K4.=| DELTAT <sob (Zs ze ie cp 


5 ZB(N+1) = ZB(N) + (KL-2.DO*K2+2.DO0*K3+K4) /6.D0 


RETURN 
END 


KRAKAKKKKKKKKRARKKKKKKKAKAKKAKRKRAKKKRAKKKKKKKAKKKKKAKAKKKKKKKKKAKARKKAKAKKAAKKA 


x SUBROUTINE STEP * 
KRRKKARKKRKRAAKKR KKK KR KKK KRK KKK KKK RRR KKK KKK KRKKKKK 


THIS SUBROUTINE PROGRESSIVELY CALCULATES THE VALUES OF ZA(N) AND 
ZB(N) BASED ON PREDICTOR-CORRECTOR CONVECTION SCHEME. 
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SUBROUTINE STEP (ZA,QVECTA,2ZB,QVECTB) 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL*8 KC,KAPPAA, KAPPAB 

COMPLEX*16 ZA(8649) ,QVECTA(8649) ,ZP,QA, QB 

COMPLEX*16 ZB(8649) ,QVECTB(8649) 

COMMON /FIRST/ ITI,ISTEPS /SECOND/ NUMPTS/THIRD/ KC,KAPPAA,KAPPAB 


DELTAT 


= 1.D0/DFLOAT(ISTEPS) 


MAKE VELOCITY EVALUATIONS FOR POINTS 1 THROUGH 4 


nom, N= 1,4 

TE = (DFLOAT(N-1) * DELTAT + TTI) 
QVECTA(N) = QA(ZA(N) ,ZB(N) ,TE) 

7 QVECTB(N) = QB(ZA(N),ZB(N),TE) 


DO 10 
NM1 = 
NM2 = 
NM3 = 
NM4 = 


=.5, NUMPTS 


N-4 


TE = (DFLOAT(N-1)*DELTAT+TTI) 


CALCULATE PREDICTED VALUE OF ZA(N) BASED ON ZB(N-1) 


i 
Z 


ZP = 


ZA 


(NM1)+(55.DO*QVECTA(NM1) -59.DO*QVECTA(NM2) 


+37. DOXQVECTA(NM3 )-9. DOXQVECTA(NM4) ) 
*DELTAT* 2.D0 * KC / 24.D0 


CALCULATE 1ST ESTIMATE OF ZA(N) BASED ON ZB(N-1) 


ie 


ZA(N) 


ZA(NM1) + (9.DO*QA(ZP,ZB(NM1),TE) + 19.DO*QVECTA(NM1) 
-5 .DOXQVECTA(NM2)+QVECTA(NM3) )*DELTAT * 2.D0 * KC / 24.D0 


CALCULATE PREDICTED VALUE OF ZB(N) BASED ON ESTIMATED ZA(N) 
ZP = ZB(NM1)+(55.DO*QVECTB(NM1) -59.DO*QVECTB (NM2) 
i +37 .DO*QVECTB(NM3)-9 .DO*QVECTB(NM4) ) 
2 ADELTAT* 2.D0 * KC / 24.D0 


CALCULATE 1ST ESTIMATE OF ZB(N) BASED ON ESTIMATED ZA(N) 


ZB(N) = ZB(NM1) + (9.DO*QB(ZA(N),ZP,TE) + 19.DO*QVECTB(NM1) 
1 -5 .DOXQVECTB(NM2)+QVECTB(NM3) )*DELTAT * 2.D0 * KC / 24.D0 


CALCULATE PREDICTED VALUE OF ZA(N) BASED ON ESTIMATED ZB(N) 
ZP = ZA(NM1)+(55.DO*QVECTA(NM1) -59.DO*QVECTA(NM2) 
i +37. DO*XQVECTA(NM3)-9.DO*QVECTA(NM4) ) 
2 *DELTAT* 2.D0 * KC / 24.D0 


CALCULATE ACTUAL ZA(N) BASED ON ESTIMATED ZB(N) 


ZA(N) = ZA(NM1) + (9.DO*QA(ZP,ZB(N),TE) + 19.DO*QVECTA(NM1 ) 
1 -5 .DOXQVECTA(NM2)+QVECTA(NM3) )*DELTAT * 2.D0 * KC / 24.D0 


CALCULATE PREDICTED VALUE OF ZB(N) BASED ON ACTUAL ZA(N) 
ZP = ZB(NM1)+(55.DO*QVECTB(NM1) -59.DO*QVECTB(NM2) - 
nt +37. DOXQVECTB(NM3)-9.DO*QVECTB(NM4) ) 
2 ADELTAT* 2.D0 * KC / 24.D0 


CALCULATE ACTUAL ZB(N) BASED ON ACTUAL ZA(N) 


ZB(N) = ZB(NM1) + (9.DO*QB(ZA(N),ZP,TE) + 19.DO*QVECTB(NM1) 
i -5 .DOXQVECTB(NM2)+QVECTB(NM3) )*DELTAT * 2.D0 * KC / 24.D0 


EVALUATE VELOCITIES AT ZA(N) AND ZB(N) 


QVECTA(N) = QA(ZA(N),ZB(N),TE) 
10 QVECTB(N) = QB(ZA(N),ZB(N),TE) 


RETURN 


END 


THIS FUNCTION SUBPORGRAM IS USED TO EVALUATE VALUES OF QA 


COMPLEX FUNCTION QA*16(ZA,ZB,TE) 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL*8 KC, KAPPAA,KAPPAB 

COMPLEX*16 ZA,ZB, DFDZ, I,AA,AB,AC 

COMMON /FIRST/ TTI,ISTEPS /THIRD/ KC,KAPPAA,KAPPAB /FOURTH/ PI 


i= 0.0, 1.D0) 
Tr = TE * 2.D0* PI 


AA = 1.D0/(ZA-1.D0/(DCONJG(ZA))) 

AB = 1.D0/(ZA-ZB) 

AC = 1.D0/(ZA-1.D0/(DCONJG(ZB) ) ) 

DFDZ =DSIN(TT)*(1.D0-1.D0/ (ZA*ZA) )-I*KAPPAA*AA-I*KAPPAS*(AB-AC) 


QA = DCONJG(DFDZ) 


RETURN 
END 


Se 


C THIS FUNCTION SUBPORGRAM DS .USHDS lose ALU ie Values Ona. 


Gi eC 4) Cl, “Oy Cle Cima 6) 0) 


* 


COMPLEX FUNCTION QB*16(ZA,ZB,TE) 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL*8 KC, KAPPAA,KAPPAB 

COMPLEX*16 ZA,ZB, DFDZ, 1I,BA,BB,BC 

COMMON /FIRST/ TTI,ISTEPS /THIRD/ KC,KAPPAA,KAPPAB /FOURTH/ PI 


I = (OepCrL. DO 
TT = TEVX 2.D0> PI 


BA = 1.D0/(ZB-ZA) 

BB = 1.D0/(ZB-1.D0/ (DCONJG(ZA) ) ) 

BC = 1.D0/(ZB-1.D0/ (DCONJG(ZB) ) ) 

DFDZ =DSIN(TT)*(1.D0-1.D0/(ZB*ZB) )+I1*KAPPAA* (BA-BB)+I*KAPPAB*BC 


QB = DCONJG(DFDZ) 


RETURN 
END 


KARARKKARKKKKKKAKRKAKKKKKKRKKAKKKAKAKKAKKKKAKKRKAKAKKKAKRKAKKRKKAKKAKKAKKAKKKKAAKKKARKRAKK 


SUBROUTINE OUTPUT * 


KRKEKKRKAKKKKKKKAKAKKAKRKKAKAKRKKAKKAKKRKRKARRARKKRKAKAKKRKKKARKAKKRKRKAKKKRKKKKKKAKRKKARKKAKRK 


THIS SUBROUTINE IS USED TO SAVE THE FINALS VALUBS SOF eZ fi) NUNG re 
FOR POSSIBLE CONTINUATION OF THE Premier 
THE VALUES ARE SAVED IN FILE “FIRRerdserooi = 
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A “Or GY ea 2 GG) 


SUBROUTINE OUTPUT (ZA,ZB) 

IMPLICIT REAL*8 (A-H, 0-Z) 

COMPLEX*16 ZA(8649) ,ZB(8649) 

COMMON /FIRST/ TTI,ISTEPS /SECOND/ NUMPTS 


DELTAT = 1.D0/FLOAT(ISTEPS) 

NN = NUMPTS - 4 

TTF = DFLOAT(NN)* DELTAT + TTI 
NUMCYC = NN / ISTEPS 


REWIND 35 
Do leeh = 1,4 
NN = NN + 1 

10 WRITE(35,*) ZA(NN) ,ZB(NN) 
WRITE(35,*) TTF 
WRITE(35,%) NUMCYC 


RETURN 
END 


RARKKKAKRKRKKRKRRRAKKRKRAKKKKKRKRKRAKRKRKRKRRRKKRKRRRKKRKRKRKRARKRKRRKAKKKKAKKRRKKKKKKKRK 


* SUBROULING GRAPH : a 


KRAKRAKRRKARKRKRKKARKRKKRKKRARKRKRKRKRKRKKRKKRKRRKKRKRKKKRKRRAKRRRRAKARARKAKKRKKKRKKKKKKRKKRK 


SUBROUTINE GRAPH (ZA,QVECTA,ZB,QVECTB) 

IMPLICIT REAL*8 (A-H, 0-Z) 

REAL*4 ZRA(8649) ,Z1A(8649) ,X(361),¥(361) ,GARRAY (1) 
REAL*4 ZRB(8649) ,ZIB(8643) 

COMPLEX*16 ZA(8649) ,QVECTA(8649) 

COMPLEX*16 ZB(8649) ,QVECTB(8649) 

COMMON /FIRST/ TTI,ISTEPS /SECOND/ NUMPTS /FOURTH/ PI 
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Cy: (OE) acGy 0s 6) 


DEFINE CYLINDER CF ey RADIUS el. 


DO 10 N = 1,361 
RAD = FLOAT(N-1) * PI /180.0 
X(N) = COS(RAD) 

10 ¥(N) = SIN(RAD) 


DEFINE REAL SIN@EE PREGYSTON NUMBEROMEOReE DOL eR GUm st ie 


DELTAT 
NPTSM4 


1.DO0/DFLOAT(ISTEPS ) 
NUMPTS - 4 


DO 20 N =1, NPTSM4 
ZRA(N) =DREAL(ZA(N)) 
ZIA(N) =DIMAG(ZA(N) ) 
ZRB(N) =DREAL(ZB(N) ) 
20 ZIB(N) =DIMAG(ZB(N)) 


GENERATE PLOTS 


CALL COMPRS 

CALL BLOWUP (1.0) 
CALL TEK618 
CALL BLOWUP (2.0) 
CALL NOBRDR 
CALL GRACE (0.0) 


VORTEN] URBCher LOT 


CALL NOCHEK 

GALL PACH (li Oye ss) 
CALL AREA2D(8.0,6.0) 
CALL YNAME('Y$',100) 
CALL XNAME('X$!,100) 
CALL YAXANG(0.0) 


CALL INTAXS 


CALL ReEKS( 2) 
CALL VILEKS (2) 
CALL GRAF(-26,1,2,-2,1,19) 


CALL CURVE(X,Y,361,0) 
GARRAY(1) = 0.08 
CALL SHADE(X,Y,361,45,GARRAY,1,0,0) 


CALL CURVE(ZRA,ZIA,NPTSM4,0) 
CALL DASH 
CALL CURVE (2ZRB,ZIB,NPISM4,0) 


CALL ENDPL(0) 


RETURN 
END 
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Figure E.lb Experimental Transverse Force Trace 
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Figure E.2) Vortex Path in Uniform Flow 
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Figure E.8 Vortex Path Z(1)=(0, 1.5084) 
T1=0.0 
(a) Cycles 1-5, (b) Cycle 6, (c) Cycles 7-15 
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Figure E.9 Vortex Paths for 100 Cycles 
Ti= 0.0 
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Figure E.10 Velocity Magnitudes for 10 Cycles 
T1i= 0.0 
(a) Z(1)=(0,.1,5083) At) eZ) aie Se) 
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Figure E.11 u-Velocity vs x-Position for 10 Cycles 
Ti=0.0 
(a) Z(1)=(0, 1.5083), (b) Z(1)=(0, 1.5084) 
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Figure E.12 v-Velocity vs y-Position for 10 Cycles 
Ti=0.0 
(a) Z(1)=(0, 1.5083), (b) Z(1)=(0, 1.5084) 
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Figure E.13 u-Velocity vs v-Velocitv for 10 Cycles 
T1=0.0 
(a) Z(1)= (0, 1.5083), (b) Z(1) = (0, 1.5084) 
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Figure E.16 Velocity Magnitudes for 10 Cycles 
Ti=0.25 
(a) Z(1)=(0, 1.6885), (b) Z(1)=(0, 1.6886) 
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Figure E.17 u-Velocity vs x-Position for 10 Cycles 
Ti= 0.25 
(a) Z(1)=(0, 1.6885), (b) Z(1) = (0, 1.6886) 


81 


“~ w wm a em ™N — oO — 
t 


PTO aaA—A 


as 


90) 


VV 


permite git MIN 


ALIOOTVIA-A 





Tee Sen 


Figure E.18 v-Velocity vs y-Position for 10 Cycles 
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Figure E.19 u-Velocity vs v-Velocity for 10 Cycles 


Ba) 2S 
(a) Z(1)=(0, 1.6885), (b) Z(1)=(0, 1.6886) 
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Figure E.22 Welocity Magnitude for 10 Cycles 
T1=0.75 
(a) Z(1)=(0, 2.2192), (6b) Z(1) =(0, 2.2193) 
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Figure E.23  u-Velocity vs x-Position for 10 Cycles 
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Figure E.24 v-Velocity vs y-Position for 10 Cycles 
T1= 0.75 
(a) Z(1)=(0, 2.2192); (6) ZO, 22298) 
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Figure E.25 u-Velocity vs v-Velocity for 10 Cycles 
Ti= 0.25 
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Figure E.26 Vortex Path Z(1)=(1.0; 1.0272) 
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Figure E.29 u-Velocity vs x-Position for 10 Cycles 
Ti=0.0 
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Figure E.30 v-Velocity vs y-Position for 10 Cycles 
T1=0.0 
(a) Z(Wie(1-0, 120272) 00S ez reo 1.0273) 
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Figure E.31 u-Velocity vs v-Velocity for 10 Cycles 
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